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Abstract 

Using an energy variational method, we calculate quasi-equilibrium config- 
urations of binary neutron stars modeled as compressible triaxial ellipsoids 
obeying a polytropic equation of state. Our energy functional includes terms 
both for the internal hydrodynamics of the stars and for the external orbital 
motion. We add the leading post-Newtonian (PN) corrections to the internal 
and gravitational energies of the stars, and adopt hybrid orbital terms which 
are fully relativistic in the test-mass limit and always accurate to PN order. 
The total energy functional is varied to find quasi-equilibrium sequences for 
both corotating and irrotational binaries in circular orbits. We examine how 
the orbital frequency at the innermost stable circular orbit depends on the 
polytropic index n and the compactness parameter GMjRc? . We find that, 
for a given GM/RcP' , the innermost stable circular orbit along an irrotational 
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sequence is about 17% larger than the innermost secularly stable circular orbit 

along the corotating sequence when n = 0.5, and 20% larger when n = 1. We 

also examine the dependence of the maximum neutron star mass on the orbital 

frequency and find that, if PN tidal effects can be neglected, the maximum 

equilibrium mass increases as the orbital separation decreases. 
04.25.Nx,04.30.Db,04.40.Dg,97.60.Jd 
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I. INTRODUCTION 



The mergers of neutron star, black hole, and neutron star-black hole binaries provide 
the most promising sources of gravitational waves for detection by laser interferometers now 
being constructed, such as LIGO and VIRGO. A variety of additional factors also motivate 
the study of binary neutron star systems: for example, such systems are known to exist with 
several having been detected in our Galaxy in the form of binary radio pulsars (eg. and 
references therein), and coalescing neutron stars may be a source of gamma-ray bursts at 
cosmological distances 0. These facts have motivated numerous theoretical studies of the 
inspiral and coalescence of compact binaries. Until recently, such calculations have involved 
either post-Newtonian (PN) point mass approximations or purely Newtonian hydrodynamic 
treatments. 

Analytic analyses of the gravitational radiation waveforms have been performed in the 
PN approximation with the binary components, whether black holes or neutron stars, being 
treated as point masses [§-0. These investigations are therefore appropriate in the early 
inspiral phase when the stars are still well separated and tidal effects can be safely neglected. 
Semi-analytic [|6|-[r3|| and numerical [|14| studies of the coalescence waveform from binary 
neutron stars also have been performed using Newtonian hydrodynamics and the quadrupole 
radiation approximation. These studies treat Newtonian effects arising from the finite size of 
the interacting fluid stars and focus on their tidal distortion, disruption, and eventual merger. 
Numerical simulations of coalescing black hole and neutron star binaries are also underway in 
full general relativity (eg. [lT5| , |l6i ). Complicated by the challenges of treating both matter and 
strong gravitational fields in (3+1) dimensional spacetime, these simulations are essential 
for obtaining definitive, quantitative results for the highly relativistic interactions which 
generate the late inspiral and final coalescence waveforms. 

There is an important distinction between binaries in Newtonian theory and those in 
general relativity: In Newtonian theory a binary can exist in true dynamic equilibrium, but 
in general relativity true equilibrium cannot exist because of the generation of gravitational 
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waves. For binaries with sufficiently large orbital separation (r >> M), the orbital decay 
time due to emission of gravitational waves is much longer than the orbital timescale, so 
that the inspiral is quasi-static through a sequence of quasi-equilibrium configurations. The 
inspiral becomes dynamic and quasi-equilibrium is destroyed with the onset of relativistic 
orbital instabilities (r ~ few M). Dynamic behavior may also commence when the purely 
Newtonian tidal instability is reached near r ~ few R (eg., For these reasons, 

both Newtonian and relativistic effects are important. 

In the relativistic simulations of Wilson, Mathews, and Marronetti |T^, the time- 



dependent equations of relativistic hydrodynamics are solved in a simplified, strong-field 
background, and the system is relaxed to find quasi-static equilibrium configurations of 
corotating neutron star binaries for a specific equation of state. These authors claim that 
neutron stars just below the maximum (isolated) mass become unstable and collapse to 
black holes when placed in close binary orbits, in contrast to the Newtonian result that 
the tidal field tends to stabilize a star (see Fig. 15 of Ref. 0; also Ref. [jl^)- Such col- 



lapse could significantly affect the inspiral and the associated gravitational waveform, and 
therefore deserves careful attention. 



Recently, Baumgarte et al. |TB| have constructed quasi-static binary equilibrium models 
in full general relativity. However, their numerical models, based on a polytropic equation 
of state, are restricted to be corotating. This work complements the numerical analysis of 
Cook [|19|, who analyzed binary black holes with zero spins. 

A fully general relativistic calculation of a binary neutron star configuration is a computa- 
tionally intensive problem. By contrast, an ellipsoidal figure of equilibrium (EFE) treatment, 
while only approximate, can find an equilibrium configuration in less than a second on a 
typical workstation. This speed affords a quick means of gaining a qualitative understanding 
of the stellar models and equilibrium sequence. These sequences can help provide a better 
understanding of the general relativistic solutions which are now being generated (eg. Ref. 
|T^). In addition, while general relativistic codes are currently limited to corotation, or 
near corotation, an EFE treatment allows straightforward construction and comparison of 



both corotating and irrotational sequences. The second order variation of the functional 
can be used to identify instabihties, e.g. the innermost stable circular orbit (ISCO) or the 
collapse of a massive neutron star in a binary to a black hole. Furthermore, both secular 
and dynamical instabilities can be identified and distinguished. 

Taniguchi and Nakamura PDI use the EFE approach to model black hole-neutron star 



binaries in the Roche approximation, using a pseudo-potential to model the strong-field 
orbital motion. Lai and Wiseman [^] insert EFE terms in the hybrid-P^N dynamical 



equations of motion, which are then integrated to identify the ISCO. The EFE approach 
has also recently been extended to first PN order in treatments of corotating binaries ||2^ . 



In this paper, we adopt the EFE energy variational method, extended to first PN order. 
Our energy functional includes terms for both the internal hydrodynamic structure as well 
as the external orbital motion. We add the leading PN correction to the internal and self- 
gravitational energy of the stars. We adopt the hybrid expressions of to PN order, to give 
the orbital terms; these expressions are exact in the test mass limit and correct to PN order 
for arbitrary masses. The energy functional is varied to find quasi-equilibrium configurations 
of both corotating and irrotational binary sequences, parameterized by orbital separation 
(or orbital angular velocity). 

In §11 we outline the basic approximations of our EFE scheme. In §111 we give the various 
energy, angular momentum, and circulation terms which comprise our functionals. The 
method of solution is discussed in §IV, and results are presented in §V. Finally, we discuss 
the implications of our results and directions for future work in §VI. In the appendices, 
we derive the PN energy contributions for an isolated spherical star obeying a polytropic 
equation of state. 

II. FUNDAMENTAL APPROXIMATIONS 

In this section, we describe the fundamental approximations and assumptions behind 
our energy variational method, which we apply to construct general Darwin-Riemann equi- 
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librium models. Further details and applications of our method can be found in ||8|-[TI . 



In steady state, an isolated, self-gravitating binary system is characterized by conserved 
global quantities such as the rest (baryonic) mass Mo and of each star and the total 
angular momentum J . The total energy of the system can be written as a functional of 
the rest mass density and velocity distributions Po(x) and v(x), respectively, even if the 
system is not in equilibrium. In principle, an equilibrium configuration can be determined 
by extremizing this energy functional with respect to all variations of Po(x) and v(x) that 
leave the conserved quantities unchanged. The fundamental approximation of our method 
is to replace the infinite number of degrees of freedom contained in Po(x) and v(x) by a 
finite number of parameters 71, 72, . . ., so that the total energy is a function only of these 
parameters: 

£; = E(7i,72,...;Mo,M;J,...). (1) 

An equilibrium configuration is then determined by extremizing the energy according to 

P = 0, ^ = 1,2,... (2) 

where the partial derivatives are taken holding the conserved quantities Mo, M^, J, . . . con- 
stant. 

An expression like equation ([^) can be written down only with the help of simplifying 
assumptions. We adopt the compressible ellipsoidal approximation to model the neutron 
stars, i.e. we assume they are triaxial ellipsoids supported in hydrostatic equilibrium by 
a polytropic equation of state [|6|,p|- [T^j23[| . Surfaces of constant density are assumed to be 
concentric self-similar ellipsoids. These approximations become exact in the limit that the 
fluid is incompressible and the gravitational potential is strictly Newtonian and truncated 
beyond the tidal term |^3[. For compressible configurations, this simplification is equivalent 



to using an ellipsoidal trial function in an energy functional and then extremizing the func- 
tional to find the equilibrium solution. Comparisons with numerical simulations show that 
the ellipsoidal model is reliable, particularly for stars obeying stiff equations of state UsHTl 
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The velocity field of the fluid is taken to be that of a Riemann-S ellipsoid, in which all 
internal velocities are linear functions of the coordinates (see eg. Ref. [§, §5, for details). The 
vorticity vector is assumed to be everywhere parallel to the orbital rotation axis. Since the 
viscosity in neutron stars is thought to be negligible, the two stars will conserve circulation 
as they inspiral so that slowly spinning neutron stars at r = oo will maintain a 

small, constant circulation. However, current general relativistic quasi-equilibrium codes 
are forced to consider corotating, or nearly corotating, sequences rather than sequences of 
constant circulation. We therefore concentrate on two interesting limiting cases: irrotational 
sequences (for non-synchronized systems with zero circulation) and corotating sequences 
(corresponding to synchronized binary systems in uniform rotation). 

Consider a binary system composed of two stars of rest mass Mq and in circular 
orbit. Throughout this paper unprimed quantities refer to the star of rest mass Mo while 



primed quantities refer to the star of mass M^. Following Ref. p3[, we denote the mass ratio 
as p = Mo/M'^. The density and pressure are related by 

P = Kp^^+^/'^\ P' = K' p'}^+^/^'\ (3) 

Note that our treatment allows the two stars to have both distinct polytropic indices {n ^ n') 
and polytropic constants {K ^ K'), although in this paper we concentrate on the case of 
two identical stars. The binary separation is denoted by r, and the principal axes of the two 
ellipsoids by ai, 02, 03, and The orientation is such that ai and a[ are measured 

along the binary axis, 02 and 03 in the direction of the orbital motion, and 03 and a'^ along 
the rotation axis. In place of the three principal axes a,, it is often convenient to introduce 
as independent variables the central rest mass density poc, and two oblateness parameters 
defined as 

Similarly we can introduce p'^^, X[, and in place of the three a[. Thus the seven indepen- 
dent variables which parameterize our models are {r, poc, Ai, A2, p'oc, A'^, Ag}. 



As is standard in EFE treatments, we assume that the density profile p{m), where m is 
the mass interior to an isodensity surface, is identical to that of a spherical polytrope with 
the same K and n, but with radius R = (aia2Ct3)^^^ [see eq. ( [B19| )]. In particular, to PN 
order, a star with semi-major axes has the same central density as a spherical star (not 
necessarily in equilibrium) with radius R in Schwarzschild coordinates, 

Mo^i GMo 1 



Po 



A7rR^\0[\ 1 Rc^ n + l 



3ai Pi {n — 3)(n + 1) 



(5) 



where we have used equations ( |A14| ) and ( |B37| ). Here 6 and ^ are the usual Lane-Emden 
variables for a polytrope (see, e.g., Ref. with 0{C,i) = 0, and where the quantities ai 
and (3i are defined in Appendix ^ and depend only on the polytropic index n (for n = 0.5 
we have ai = —3.39307 and (3i = 3.63708, while for n = 1.0 we have ai = —2.77960 and 
Pi = 2.67168). In the Newtonian limit (c — > oo), equation (|) reduces to the expression 
used in the Newtonian treatment of |]8|-[TT|. Although equation (|^) originates from a PN 
expansion, our code treats this relation as exact. 



III. ENERGY, ANGULAR MOMENTUM, AND CIRCULATION TERMS 

A. Newtonian Treatment 

We now review the energy, angular momentum, and circulation terms in the Newtonian 
EFE approximation, which has been covered in detail in [P-p!l|. When similar expressions 
can be written for the two stars, we give only the one corresponding to Mo, the other being 
obtained simply by replacing unprimed with primed quantities. 

The usual Newtonian orbital contributions to the energy and angular momentum are 

1 2 2 gmqm:, , 

En orb = T^P^l r , JNorb = pr il, (6) 

2 r 

where fj, = MoM'^/ {Mo + M^) is the reduced mass and VL = v/r is the orbital angular 
velocity. The gravitational tidal (quadrupole) interaction energy Wt between the two stars 
is (see Appendix B of Ref. §) 
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W, = - I',, - /^3) - ^(2^1 - h2 - hs). (7) 



Here we have defined 



I^^ = M^IiMm^ Vo"c'/' (no summation) , (8) 



where 



hi = yI' ^2 = 7^, h = A1A2, (9) 



with analogous relations for Jj^ and h'j. We will refer to key relations in Ref. |^ by their 
equation number preceded by an "I"; for example, equations (||) and @ are also equa- 
tions (I.B8) and (I.B9), respectively. The dimensionless structure coefficient ^3 appearing 
in equation (||) is defined by [cf. eq. (1.3.17)] 

where 

so that K„ = 1 for n = 0. Values of and k„ for various n are tabulated in Table I of 
The internal energy of the star with mass Mq is given simply by 

U = hKpli'^Mo. (12) 
The self-gravitational energy can be written [cf. eq. (1.4.6)] 

W = -hGM'J'pli'f, (13) 
where the dimensionless ratio / is given by 

^ ~ = 2(a,a2a3)V3 ' (14) 



so that / = 1 for a spherical star. The index symbols Ai are defined as in Ref. [^] (§17), 



y4j = 010203 / — — ^ -, with A'^ = {al + u){al + u){al + u), (15) 

JO /\[a^ + u) 

and are functions only of Ai and A2. In equations (|l^) and (p!3D , ki and k2 are dimensionless 



structure constants depending only on the polytropic index n, defined as 

Values of ki and k2 are tabulated in Ref. §]. 

We permit the star to be differentially rotating with uniform vorticity (e^ as measured 
in the corotating frame of the binary: 

C ^ (V X u) . ea = -^i±^A. (17) 

0102 

Here the velocity field u in the corotating frame is given by 

u = QiX2ex + Q2Xie2 (18) 

with 

Qi = -^^C = +-A, 

af + 02 02 

Q2 = +^^C = --A. (19) 
af + 02 oi 

The quantity A is the angular velocity of the internal fluid motions in the corotating frame. 
Here ei is along the binary axis, directed from Mo to M^, 62 is in the direction of the orbital 
velocity, and 63 is perpendicular to the orbital plane. The velocity field in the inertial frame 
is given by 

u(°) = u + O X X, (20) 

and the vorticity in this frame is 

C(o) = (V X uW) ■ es = (2 + fn)n, (21) 

where 
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Corotation (synchronization) corresponds to /ij = C = A = 0, while an irrotational velocity 
field is obtained when = —2. Note that the geometric outer shapes of the two stars 
always rotate at the orbital angular velocity Q in the inertial frame, regardless of the internal 
velocity field. 

The Newtonian spin kinetic energy Tg (i.e., the kinetic energy in internal fiuid motions) 
in the inertial frame can be obtained from equations (p!8|)-(pOD. Assuming a Lane-Emden 
trial density function, one finds [cf. eq. (1.5.6)] 

T, = h{A^ + n^)-2 (/n/22)'/' Afi, (23) 

where / = In + 122 is the moment of inertia. Similarly, the Newtonian spin angular momen- 
tum Jg can be written [cf. eq. (1.5.5)] 

Js = m-2{Iul22f^A. (24) 
As given by P-pJ]], the purely Newtonian total energy of the system, not necessary in 



equilibrium, is therefore 

En = Em orb + U + U' + W + W' + Ts + T', + Wt. (25) 
Similarly, the total angular momentum to Newtonian order is 

Jn = JNorb + Js + J's- (26) 

Another important quantity, conserved in the absence of viscosity, is the fiuid circulation 
C along the equator of the star (see §5.1 of Ref. P]). We define 



C 



-^K^Mo) C = IA-2 (/n/22)'/' ^. (27) 
on / 



The quantity C has the dimensions of angular momentum but is proportional to the conserved 
circulation C. For convenience, we typically refer to C itself as the circulation. 
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B. Post-Newtonian Treatment 



In this paper, we augment the Newtonian terms of § 111 A with PN contributions (a) to the 
orbital energy and angular momentum, (b) to the internal and self-gravitational energy of a 
star, and (c) to the coupling energy and angular momentum between the internal structure 
and the orbital motion. 

We replace of the Newtonian contributions to orbital energy and angular momentum [eq. 
(^)] with the hybrid-PN expressions of 0] (with r = 0, valid for circular orbits): 



E 



H 



Et + Er,, Jh = Jt + Jrj 



where the Schwarzschild, or test-mass, contributions are (with G = c 
equations only) 



Et = fi 



r — m 



r + m 



1/2 



(r + mY 



r — m 



-1/2 



f 


(r + m)^ 


1/2 - 




r^(r — m) 





(r + m)" 



- 1 



-1/2^ 



r^(r — m) 

and the PN corrections unaccounted for by the Schwarzschild terms are 

'9v^ lGm\ 



(28) 

1 for these two 

(29) 
(30) 



En 



8c^ 



2 rc^ 
Gm\ 



rc 



(31) 
(32) 



Here v = Qr is the orbital velocity, the total rest mass m = M^+M^, and the mass parameter 
7] = fi/m. In the notation of Ref. IQ, m also includes contributions from the self energy of 
the stars and is not simply the total rest mass; however, we account for these contributions 
below [see eqs. (^) and (0)] as a coupling between the internal structure and the orbital 
motion. With our definitions, the quantities m and rj remain strictly constant along the 
sequence. The radius r is the so-called harmonic or de Bonder coordinate. When expanded 
in Gm/rc^, the above expressions for E^ and Jh give the Newtonian terms [eq. (^)] to 
lowest order and the correct orbital PN terms to first order. 
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The PN contributions to the internal and self-gravitational energy of a spherical poly- 
trope are derived in Appendix ^: 

UpN = -h^Kpt^'M^/^ (33) 
Wp^ = -h^Mll\-i\ (34) 



where the coefficients l\ and 1% are defined by equations ( [B42|) and ( [B52|) , respectively, and 



depend only on the polytropic index n. Values for l\ and are listed in Table |. We 
also include the analogous expressions for the PN energy contributions f/p^ and Wp^ in the 
second star. For isolated stars, variation of the energy function E{poc) = U +W +Upn + Wpn 
does a good job of approximating the exact equilibrium sequence obtained by integrating 
the Oppenheimer-Volkoff (OV) equation, especially for soft equations of state (see Figure p. 

We neglect PN corrections due to asphericity; for small perturbations from spherical, 
equations ( p3D and (|3^ are valid to lowest order. To account completely for all PN cor- 
rections to the internal and self-gravitational energy, equations (^) and ( p4[ ) would need 
to be multiplied by a form factor analogous to /(Ai, A2) in equation ([T^). We also neglect 
P^N and higher order corrections to the internal and self-gravitational energy of the neutron 
stars. 

The final PN correction terms we include involve the coupling between the orbital motion 



and the internal structure. Effacement of internal structure (see Ref. and references 
therein) requires that we substitute the effective masses M^ff = Mo + Egeif/c^ and M'^j-j- = 
M'o + E'^eifl Tilace of the rest masses Mo and in the Newtonian equation (P) if we 
wish to maintain PN order accuracy. Here the self energy of the mass Mo neutron star is 
Eseif = W + U + UpN + WpN- We have chosen to exclude the spin kinetic energy Tg from 
the self energy E^^if in order to ensure that our binary configurations satisfy the identity 
dE = QdJ exactly (see §fV|). We speculate that additional higher order terms which we are 
neglecting would restore this identity when we include Tg. Furthermore, Tg is typically small, 
especially for irrotational sequences, so that we are neglecting a term which is smaller than 
our order of accuracy. Expanding to PN order, the resulting corrections to the Newtonian 
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point-mass energy and angular momentum are 



Eseif (lMl^_GK 



) 



(35) 



\2 IV? r 




(36) 



with analogous expressions E'qj and J'qj due to the corrections arising from the self energy 
of the neutron star of mass M'^. The subscript 01 serves as a reminder that these terms 
couple the orbital motion to the internal structure of the neutron stars. Although including 
f/pAT and Wpj^r as part of Eg^if in equations ( P3[ ) and (^) amounts to keeping terms which 
are formally of P^N order, doing so keeps our energy functional accurate even for stars at 
r = oo but with t> > 0; excluding Upn and Wpn from Egeif would lead to an energy functional 
which yields manifestly incorrect results for isolated stars (in particular, the stability of a 
star would depend on the velocity of the observer). 

We neglect PN effects in the internal motion which would appear as corrections to the 
spin kinetic energy Tg, the spin angular momentum Js, and the circulation C. For corotating 
sequences, these corrections are on the order of R/r times the PN orbital terms, so that we 
are justified in neglecting PN internal motion effects all the way down to roughly the ISCO 
(r ~ few R). For irrotational sequences, PN internal motion is especially small, since the 
fluid motion in the inertial frame is always predominately set by the orbital motion. We are 
also neglecting all PN tidal terms (but see § |V1 C|) . 

The total energy and angular momentum functionals used in our treatment are therefore 



E = Eh + U + U' + W + W' + Ts + T^ + Wt + UpN + U'pj^ + Wpn + IVp^v + Eqi + E'^ 



Since we are approximating the internal fluid motion as Newtonian, we continue to use 
equation (|27|) for the circulation. 



(37) 



and 



J — Jh + Js + J's + Joi + J'oi- 



(38) 
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IV. FINDING EQUILIBRIUM SEQUENCES 



In order to find an equilibrium solution, we must extremize the energy function E while 
holding all the quantities Mq, M^, C, C, and J constant. We are therefore identifying those 
configurations for which the energy is unchanged to first order by any perturbation which 
conserves rest mass, circulation, and total angular momentum. An equilibrium sequence is 
constructed by repeating this extremization at various values of the orbital separation r. In 
the case of corotation, the stars are members of the Riemann-S family with vorticity C = 0, 
and although the perturbations considered conserve circulation, the circulation does indeed 
vary along the equilibrium sequence itself 



Unfortunately, the various energy terms in §|T| do not depend explicitly on J, C, 
and C. Instead, the energy has been written in terms of the variables from the set 
{r, poc5 Ai, A2, Poc5 -^'i' -^25 ^5 ^5 ^'l' thcsc varlaHcs are not independent. The Lagrange 
multiplier method allows us straightforwardly to minimize the total energy E{ai) subject 
to the constraints that the total angular momentum J(crj) and the circulations C((jj) and 
C'((Ti) are constant. The equations which determine an equilibrium configuration are 

dE dJ dC ^ dC 

■J— + Ojr- + ^CTT- + ^c'T^= 0, 39 
acTj acTj oai ocTi 

where ^j, ^c, and ^c' are the Lagrange multipliers. Here the partial differentiation with 
respect to cTj is done with all the remaining variables aj {j 7^ i) being held constant. Although 
equations (^) are accurate only to the order of our expansion of E, J, C, and C, our 
numerical code treats them as exact. We expanded and restructured a version of the code 
developed in Ref. [|ll|] to solve these coupled equations numerically. For a given value of 
r, we use a Newton-Raphson scheme to solve the 10 equations (|39D and find the variables 
Poc, Ai, A2, p'oc, A'l, A2, ^, O) ^c, and ^c'- In practice, our code uses Cramer's rule on the ctj G 
{Q, A, A'} equations of (|39D to find the Lagrange multipliers: 
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dE dC dC 

an an an 

as ac ac 

dA dA dA 

dE ac ac 

aA' aA' aA' 



aj dE ac 

an an an 

dj dE ac 

aA aA aA 

dj dE ac 

dA' dA' dA' 



dJ 


ac 


dE 


an 


an 


an 


aj 


ac 


dE 


aA 


aA 


aA 


dJ 


ac 


dE 


aA' 


dA' 


dA' 



(40) 



where the determinant 



aj ac ac 

an an an 

aj ac ac 

aA aA aA 

aj ac ac 

dA' dA' dA' 



(41) 



Using equations (^0]), we substitute for the Lagrange muhiphers in the remaining 7 equations 
of ( p9D to find the orbital velocity as well as the shape and structure of the neutron stars. 

The Lagrange multipliers can even be determined analytically, since expressions for the 
various derivatives appearing in equations (1^) and (^Tj) can be derived from equations (pTj), 
(p8|), and (pTj). The algebra is simplified somewhat by dC/dA' = and dC /dA = 0. We find 
O = (,c = and ^c' = as mathematical identities, for any r, poc, Ai, A2, p'^c, ^i, ^2 
and not just along equilibrium sequences. This provides a convenient check of a large portion 
of the code, since the Lagrange variables are determined by evaluating the determinants in 
equations (l40| ) and ( ^if ) directly at every iteration in the Newton-Raphson scheme. Further- 
more, equation ( ^9|) then implies 



dE - VLdJ - AdC - A'dC = 0, 



(42) 



so that dE = QdJ along both corotating (A = A' = 0) and constant circulation [dC = 
dC = 0) sequences. Consequently, the minima of J and E always coincide precisely on our 
equilibrium sequences. 



V. RESULTS 

We concentrate on binaries composed of identical stars, each with polytropic index n = 
0.5 or 1.0, and with a compactness parameter GM/Rc^ = 0.125,0.2, or 0.25 in isolation. 
These values of n and GM/Rc^ cover an interesting region of parameter space relevant 
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to neutron stars. For the AV14+UVII equation of state of Ref. [^], stars in the mass 
range 0.7 Mq ^ M < I.TMq have compactness parameters 0.1 ^ GM/Rc^ < 0.25 with 
corresponding effective polytropic indices 1 ^ n ^ 0.5 (see Table 3 of Ref. [|10[)- For an 



isolated neutron star with M = 1.4Mq, the nine equations of state in Ref. which allow 



this mass all yield a compactness parameter in the range 0.13 ^ GM/Rc^ ^ 0.28. 

Some physical quantities describing our neutron stars in isolation are listed in Table J|. 
These quantities can be converted to physical units once a value of the polytropic constant 
K is chosen. For instance, for n = 0.5 and K = 1.35 x 10~^° cm^/g^/sec^, the maximum 
gravitational mass M""''^ = 2.2OM0, so that M = 0.83, 1.39, and 1.74Mq and R = 9.8, 10.3, 
and 10.3 km for GM/Rc^ = 0.125, 0.2, and 0.25, respectively. Forn = 1 and K = 1.20 x 10^ 
cmVg/sec^, the maximum mass M™"^ = l.GOOM© so that M = 1.15,1.52, and 1.599M0 
and R = 13.6,11.2, and 9.4 km for GM/Rc^ = 0.125,0.2, and 0.25, respectively. Any of 
our models can be scaled to a particular mass (say I.4M0) by an appropriate choice of the 
polytropic constant K. 

Integration of the OV equation reveals that the true maximum masses of n = 0.5 and 
1.0 sequences occur at log Kpoc/c^ = —0.082 and —0.497, respectively, at which point 
GM/Rc^ = 0.316 and GM/Rc^ = 0.214. We note, however, that our n = 1.0 isolated 
stars with GM/Rc^ = 0.25 are indeed stable in our PN-accurate model, since our energy 
functional [eq. ( P44|) ] gives GM/Rc^ = 0.2556 at the maximum isolated mass (see Table 



OD- 

When expressing our results, the orbital angular velocity, angular momentum, energy, 
and circulation are listed in terms of the dimensionless ratios 

(7rGp„)V2' (GM3i?,)i/2' {GM^JRS {GMlR^y/^' ^ ^ 

where po = Mo/(47ri?^/3). Here Ro is defined to be the radius that an isolated polytrope, 
with the same Mq, K, and n, would have in Newtonian theory (ie. without the Upn and 
WpN corrections): 

'{n+l)K 

AttG I \^<i\0'i 



Ro = ^1 



n/(3-n) / X {l-n)/(3-n) 



2,..,/ • (44) 
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The radius Ro is a convenient unit of length because it is constant along equilibrium se- 
quences. 

The orbital separation r is measured in harmonic coordinates while the semi-major axes 
Oj are measured in Schwarzschild coordinates, which is permissible because we have two 
distinct small parameters in our PN expansion, Gm/rc^ and GM/Rc^. We determine when 
the neutron stars touch by the PN accurate relation r + GMt/(? = ai + a'^, where Mt = 
2Mo + E/c^ is the total mass energy of the system. Since the radial coordinate r is not a 
gauge invariant quantity, we choose to make our plots as a function of the orbital angular 
frequency / = Q/2'k (= 1/2 the gravitational wave frequency). The non-dimensionalized 
orbital velocity Q can be converted to a frequency / in Hz via the relation 

M,,/^ 1788H.(5^) (45) 

where M1.4 = M/1.4Mq. For given values of n and GM/Rc^, the corresponding values of 
GMo/ Ro(? and M /Mq can be determined from Table 0. For example, an n = 1 model with 
GM/Rc"^ = 0.2 corresponds to GMa/RoC^ = 0.1467 and M/M^ = 0.9083. While and 
are constant along our equilibrium sequences, the ratio GM/ Rc^ does indeed vary; however, 
whenever we specify a particular value of GM/ Rc^ we are referring to that value at infinity 
(ie. for isolated neutron stars). 



A. Irrotational Sequences 

Figures ^ and |^ present irrotational (C = 0) sequences in the equal mass case with 
n = 0.5 and n = 1, respectively, showing the total energy E, relative to its value E^o at 
infinity, for GM/Rc^ =0.125, 0.2 and 0.25 (thick solid curves). The energy E^o can be 
calculated from equation ( |B44| ) for each star. Also shown in this figure are the point-mass 



Newtonian equilibrium energy (thin dashed curve), as well as the results of a Newtonian 
EFE treatment as in P|-pTJ (thick dashed line). Figure ^ shows the dependence of GM/ Rc^ 
on the orbital frequency fisco at the ISCO for n = 0.5 and n = 1 irrotational sequences. 
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For each n, the curves are terminated at the value GM/ Rc^ corresponding to the maximum 
mass. 

Table |T1| presents physical quantities along selected irrotational sequences for n = 0.5 
and 1.0, and various values of GM/Rc^. For each sequence, the value of r/ai marked by a 
single asterisk (*) identifies the equilibrium configuration for which E (and J) is a minimum, 
ie. the ISCO. As expected, the size of the neutron star R/Ro decreases as GM/Rc^ increases. 
Our GM/ Ro(? = sequences agree with the Newtonian results of Ref . 0] (see their Table 

1). 

Tables |V| and show the sizes of the individual energy and angular momentum terms 



appearing in equations (|37D and (^) for the specific sequence GM/Rc^ = 0.2, n = 1. Note 
that some of the relativistic corrections are included in and Jh'- for instance, at the 
ISCO, Eh is 1.7 times smaller than what the Newtonian expression [eq.(^] would give. The 
importance of the PN terms increases as GM/ Rc^ increases, and the stars behave more like 
point masses. 

B. Corotational Sequences 

Fully relativistic calculations of binary neutron star systems (eg. |ll6| , [T8| ) provide a nice 



test of our PN models. However, such calculations are currently limited to corotating, or 
nearly corotating, stars. Figures ^ and ^ present corotating (A = A' = 0) sequences with 
n = 0.5 and n = 1, respectively, showing the total energy E, relative to its value E^o at 
infinity, for GM/Rc^ = 0.125,0.2, and 0.25 (thick solid curves). 

Tables presents physical quantities along selected corotating sequences. For each 
sequence, the value of r/ai marked by a single asterisk (*) identifies the equilibrium con- 
figuration for which E (and J) is a minimum, the innermost secularly stable circular orbit, 
while the value marked by a double asterisk (**) identifies the innermost dynamically stable 
circular orbit. The dynamical ISCO occurs when = dE/dr at constant circulation along 
an equilibrium sequence (see Ref. We numerically evaluated dE/dr simply by finding 
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equilibrium configurations at two slightly different orbital separations and then finite dif- 
ferencing. As in our irrotational sequences, the overall size of the neutron star R/Rq does 
indeed decrease as GM/Rc^ is increased. For corotation, the spin of each star causes a 
rotational bulge such that 02 > 03. The circulation C grows in magnitude as the radial 
separation of the neutron stars decreases. 



C. Innermost Stable Circular Orbit 



Table |VII| presents the total angular momentum J, the total mass energy Mt = 2Mo + 
E/c^, a^jaxi O3/01 and the orbital frequency /, all at the ISCO of selected irrotational 
sequences. For any given n, the orbital frequency / can be converted into units of Hz if 
either the gravitational mass M or the equation of state (ie. the polytropic constant 
is specified. It is therefore useful to express the ISCO frequency both as M\Afisco [Hz] 
(when M is known) and M^^^ fisco [Hz] (when the equation of state is known), where 
M1.4 = M/{1AMq) and M]^^'' = M"'"'' / (I.GMq). Here M"""^ is the maximum gravitational 
mass in isolation for a given equation of state. The conversion factor between Mi^fjsco 
and M^Q^ fisco depends on the quantity M/M™"^, which can be determined from the last 
two columns of Table |T|. 

Table [VIII| lists the orbital frequency ffsco innermost secularly stable orbit as well 

as the frequency ffsco innermost dynamically stable orbit. For the six combinations 

of n and GM/Rc^ listed in this table, the frequency ffsco differs by less than 4% from the 
fisco values of the corresponding irrotational sequences. That the spin has only a small 
effect on the orbital dynamics is consistent with the small values of Ts/|l^| (typically on the 



order of 0.01 near the ISCO) listed in Table Wl 



D. Maximum mass 



An important issue is whether the maximum mass of each neutron star increases or 
decreases as the binary components are brought together. Figure |^ shows the rest mass 
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Mo, normalized to the maximum rest mass M'^'™^^ at infinite separation, as a function 
of tlie central density parameter Kp\l^ j(? for an n = 1 corotating binary system and at 
three different values of the orbital angular frequency: / = (solid curve) corresponds to 
spherical isolated stars, M"^^^ f = 434 Hz (long dashed curve) corresponds to the innermost 
dynamically stable orbit for a star with Mo/M^""'"^^^ = 0.938 {GM/Rc^ = 0.2 at infinity), 
and M^g^f = 463 Hz (short dashed curve) corresponds to the innermost dynamically stable 
orbit for a star with Mo/M^''^'''' = 1 {GM/Rc^ = 0.2556 at infinity). As they inspiral, stars 
which begin on the stable side {dMo/dpoc > 0) of the / = curve move along horizontal 
lines (since Mo is conserved) to the left on this plot. 

Figure ^ shows the maximum equilibrium rest mass M™*^^ as a function of the orbital 
angular frequency for corotating (dashed curves) and irrotational (solid curves) sequences. 
For example, the maxima of the three curves in Fig. ^ gives three data points for the n = 1.0 
corotating curve in Figure The curves are terminated at the values of / for which stars 
of rest mass = ]Vf^'™^^ acquire a dynamically unstable orbit. We see that the maximum 
equilibrium mass increases as the orbit decays, regardless of the sequence considered. We 
conclude that if the maximum equilibrium mass does ever decrease in a binary system, it 
can do so only through PN tidal terms, or higher order PN terms, which we are neglecting 
in this paper (see ^VI C|) . Note that we have not determined from our equilibrium analysis, 
which of our equilibria are stable. The exception is the / = (infinite separation) sequence, 
for which the turning point on the M versus Poc plot marks the onset of radial instability to 
collapse. 

VI. DISCUSSION 

A. Summary of Results 

We have extended the work of previous studies in a number of ways. We have improved 
upon purely Newtonian EFE treatments by including PN effects, and improved upon point- 
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mass PN treatments by including both Newtonian and PN finite-size effects. The main 
message of Appendix B is that equation ( |B44| ) is the PN accurate energy functional of 
an isolated non-rotating star of arbitrary polytropic index n. Nowhere in the derivation 
of equation ( [B44| ) is equilibrium assumed, but in equilibrium the trial density function we 



implement agrees to PN order with the density profile obtained by integrating the OV 
equation. The PN coefficients li and I2 appearing in equation ( |B44D are listed for various 
polytropic indices n in Table |. Figure |l] demonstrates that using equation ( [B44|) yields 
reliable results (compare the long dashed curve to the OV results given by the solid curve), 
especially for soft equations of state (which have a maximum mass at small poc)- One 
particularly nice feature of our energy functional is that it correctly exhibits a maximum 
mass even for stiff equations of state. Various physical quantities for our isolated spherical 
models are listed in Table ||. 

Figures ^ and H, present the equilibrium energy £" as a function of the orbital frequency 
/ along irrotational sequences with various values of the polytropic index n and compactness 
parameters GM/Rc^. Note that for a given M and GM/Rc^, n = 1 polytropes are more 
centrally concentrated than n = 0.5 polytropes, so that the former behave more like point 
masses (ie. Newtonian effects which destabilize the orbit are less important for n = 1) and 
have a slightly larger Mi^fisco- This same tendency can be seen by comparing the n = 0.5 
and n = 1.0 corotating sequences presented in Figures ^ and |^, respectively. In these plots 
the minimum of the equilibrium energy E marks the onset of an orbital secular instability. 
The dynamical instability sets in at higher orbital frequency, as can be seen from Table 

or [vmi . 



Tables |ITT| and ^ present various quantities along selected sequences for comparisons 
with future studies. The Tables [V| and [V| present the relatives sizes of the terms which 
make up our energy and angular momentum functionals; these tables allow us to explicitly 
examine the relative importance of the various effects which determine the orbital dynamics 
and the interior structure. Comparing the magnitudes of Wt and Eqi, it is clear that 



both tidal and relativistic terms play an important role. Tables |VI1| and |VIII 
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summarize 



some of the important quantities at the ISCO of our irrotational and corotating sequences, 
respectively. Note that for fixed n, Mi^^fisco always increases with GM/Rc^, since more 
compact stars behave more like point masses. 

Although real neutron star sequences are probably close to irrotational, fully relativistic 
calculations of binary systems are currently limited to corotating sequences. By looking 
for turning points in the equilibrium energy, it is straightforward to determine the secular 
ISCO frequency ffsco ™- ^^^^^ corotating sequences. It is often assumed that the actual 
ISCO frequency fisco of an irrotational sequence does not differ drastically from the fre- 
quency ffsco determined from corotating calculations. Our results allow us to quantify this 
difference. For n = 0.5 we find that the irrotational fisco is approximately 17% larger than 
the corresponding ff^coi while for n = 1 the difference is approximately 20% (compare 
entries in Tables [VI I| and |VIII|) . For any polytropic index n, this difference depends only 
very weakly on GM/Rc^. 

If multiple binary neutron star coalescence events are observed by instruments such as 
LIGO, then tables like Table |Vli] , or figures like Figure ^ will give information regarding the 
equation of state. For instance, consider a simple world in which binary components always 
had identical masses, although allow for varying neutron star mass from one binary to the 
next. Suppose further that all neutron stars obeyed the same (but unknown) polytropic 
equation of state. A single coalescence observation would yield a value of Mi ^fjsco: since 
the gravitational mass M is encoded in the early inspiral waveform and since the ISCO 
frequency fisco can be determined by the onset of the orbital plunge. However, the equation 
of state and compactness parameter GM/ Rc^ could not be uniquely determined from a single 
observation, since we do not yet know which curve (if either) in Figure ^ applies. For each 
value of n being considered, the corresponding maximum gravitational mass M""*^ can be 
immediately determined from the observed value of fisco and from the M^q'^ fisco column 
in Table |V14 An additional observation of a coalescence event with a different fisco would 
yield another value of M™'"^ for each n. The polytropic index which consistently gives the 
same M"^°-^ regardless of the observed fisco must then represent the true equation of state. 
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Tables similar to Table [VI1| will therefore be helpful in testing realistic equations of state 
once coalescence observations are available. 

There is only one equation of state for neutron stars, even if knowledge of it is still 
uncertain. For a given equation of state, one might ask how fisco depends on the mass of 
the binary components. For polytropes, fixing n and K also fixes the maximum gravitational 



mass M™"^^ in isolation. From the M'^q^ fjsco column of Table |VI1| , we see that for n = 
0.5 the general trend of fisco is for it to decrease as the stellar mass M (or GM/Rc^) 
increases, while for n = 1 the frequency fisco increases with M. We therefore conclude 
that the relation between fisco and M depends strongly on the equation of state. The 
reason for this dependence can be traced to two competing effects. Namely, as we consider 
more compact (larger M) stars, PN orbital destabilization effects become stronger, while 
Newtonian destabilization effects become weaker. 

For the cases listed in Table [VI H , the total angular momentum J is 0.96 to 1.10 times 
the total energy squared M"^ = (2Mo + E)'^ (in gravitational units, G = c = 1) at the ISCO. 
In those cases with J > M^, angular momentum would need to be either radiated away or 
carried away by ejected matter before a black hole can form (see Ref. p3] for a discussion). 



We find J > for all n = 0.5 sequences with GM/Rc^ < 0.21, while for n = 1 this 
criterion is satisfied when GM/Rc^ ^ 0.18. 

Figure ^ presents the equilibrium rest mass Mq as function of the central density for 
various orbital frequencies / for an n = 1 corotating sequence. We see that, for the terms 
we have included in our energy functional, the maximum mass in a close binary is slightly 
larger than in isolation. Figure | summarizes our maximum mass results for both irrotational 
and corotating sequences, and for n = 0.5 and 1.0. The irrotational sequences maintain a 
lower maximum mass than their corresponding corotating sequences, because corotating 
stars (a) have more rotational support and (b) are more ellipsoidal so that the Newtonian 
tidal field has more of a stabilizing effect. The relative increases in M™"^ is more gradual 
for n = 0.5 than n = 1.0, because n = 0.5 stars near the maximum mass are much more 
centrally condensed than n = 1.0 stars near their maximum mass: from Table f| we see 
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that both the central density parameter Kpod (? and the compactness parameter GMj R(? 
are much larger for the maximum mass n = 0.5 star (log Qc = 0.12, GM/Rc^ = 0.39) than 
the n = 1 star (log Qc = —0.39, GM/Rc^ = 0.26). Consequently, the n = 0.5 stars near 
the maximum mass have smaller moments of inertia and behave more like point masses, so 
that stabilizing hydrodynamic effects play a less significant role. Our results suggest that 
the maximum equilibrium mass increases as the orbit decays, regardless of the sequence 
considered. However, if the maximum equilibrium mass does decrease in a binary system, 
it must be a result of the PN tidal terms, or of higher order PN corrections, that we are 



neglecting (see § [VI Q) 



B. Comparison with Other Work 



Recently, Taniguchi and Nakamura |^ have applied the ellipsoidal approximation to a 
neutron star-black hole binary. They also compute irrotational and corotating sequences, 
but only for n = polytropes. In place of the hybrid-PN orbital energy contribution 
Eff, they use a generalized pseudo-Newtonian potential constructed to fit the ISCOs of 
the hybrid P^N treatment of Kidder, Will, and Wiseman [Q. We improve upon the work 
of Taniguchi and Nakamura by including the relativistic corrections Up^ and Wp^ to the 
internal energy of the neutron star, and by calculating sequences of fixed Mo/Ro, instead of 
fixed M/R which varies as two neutron stars inspiral. Despite these differences our results 
do agree qualitatively: for a given equation of state more compact stars have larger ISCO 
frequencies, and for a fixed compactness parameter the ISCO along irrotational sequences 
is larger than the secular ISCO along corotating sequences. 



Taniguchi and Shibata |2^ have recently presented PN EFE analyses of binary neutron 
star systems which is rather similar to our own. The main differences are that (1) they derive 
equilibrium conditions from a tensor virial method, whereas we use an energy variational 
method, (2) they include all terms, up to quadrupole order, which are formally PN order, 
while we have neglected PN tidal terms and PN corrections to the internal fluid motion, (3) 
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they consider only corotating sequences, while we also treat the more realistic irrotational 
sequences, and (4) they implement a simple Lane-Emden trial density function, as opposed 
to our PN-accurate density function [see eq. |B19|| . Our treatment does somewhat better 



than Taniguchi and Shibata's at matching the numerical, PN-accurate results of Shibata 
7^. For instance, for n = 0.5 and GMo/RqC^ = 0.02, Shibata finds secular ISCO values 



of E, J, and equal to —1.236, 1.457, and 0.31, respectively; compare this to our values 
—1.246, 1.458, and 0.316, and to Shibata and Taniguchi's values —1.269, 1.443, and 0.327. 
We find a similar level of agreement with Shibata along other sequences: for n = 1.0 and 
GMo/RqC^ = 1/60, Shibata finds secular ISCO values of E, J, and VL equal to —1.146, 1.386, 
and 0.35, while our treatment gives —1.155, 1.384, and 0.364. 

In order to examine the ISCO in irrotational sequences, Lai and Wiseman have 



combined the Newtonian ellipsoidal equations of motion |[T^ with the hybrid P^N point 
mass equations of motion of Kidder, Will, and Wiseman. This treatment does not include 
the PN coupling between the orbital motion and internal structure {Eqi and Jqi in this 
paper), nor does it include PN corrections to the internal and self-gravitational energy (f/pAr 
and Wp]s[ in this paper). Nevertheless, our results are in excellent agreement with those 
of Lai and Wiseman (see the comparison in Table |Vli| ): for the three values of GM/Ro(? 
for which a direct comparison is possible, we find agreement to better than 4% in fisco 
and to better than about 2% in the ratios 02/01 and a^/ai. The comparisons of 02/01 and 
of 03/01 should be viewed with some caution because of the freedom to measure distances 
with various coordinate systems in general relativity. In the Newtonian limit our results 
agree precisely with the Newtonian results of Lai and Wiseman [see their equations (8) and 
(9)], as they must since both treatments are extensions of the same Newtonian ellipsoidal 
approximation. 

In the point mass limit, the results of Lai and Wiseman are accurate to P^N order, 
whereas our results are accurate only to first PN order. We have not implemented the 
point mass P^N orbital energy and angular momentum corrections in our energy functional, 
because doing so causes dE = QdJ to hold to P^N order but not exactly along our sequences, 
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and because our functional would be neglecting other P^N terms which are larger. However, 
we have explicitly tested the importance of the point mass P^N corrections by adding them 
to our energy and angular momentum functions. We find that for GM/Rc^ = 0.25, the 
ISCO frequency for our n = 0.5 and n = 1 irrotational sequences is decreased by less than 
4%, while for GM/Rc^ = 0.125 and 0.2 this frequency is increased by about 2% or less. 

Baumgarte et al. |jl8| have recently calculated in full general relativity quasi-equilibrium 
corotating sequences of binary neutron stars obeying a polytropic equation of state. Their 
computations are the most reliable to date for the neutron star structure, maximum equi- 
librium mass, and determination of the secular ISCO. Over the range of orbital separations 
they consider, they find the maximum equilibrium mass to increase slightly as the orbit 
decays, in agreement with the results presented in Figures |^ and |^. This agreement is reas- 
suring, since these relativistic calculations essentially represent the "true" solution. Figure |^ 
compares our binding energy results with those of Ref. [|1^ for an n = 1 corotating sequence. 



The maximum equilibrium rest mass results presented in Figures |^ and ^ are also in 
qualitative agreement with the results of Lai [|T^: the maximum mass of a neutron star 
seems to increase as the orbit decays. The main differences between our energy functional 
and that of Lai are that (a) Lai's adopted equation of state is based on a degenerate, non- 
relativistic Fermi gas of nucleons (n = L5 polytrope), while we have concentrated on n = 0.5 
and n = 1 polytropes, (b) Lai has included a correction AEi^t as in Ref. |0 accounting for 



the tendency of neutrons to become mildly relativistic, and (c) we use equation ( |B44| ) for 
the self energy of a polytrope, whereas Lai uses equation ( |B47| ) . Our maximum mass results 
are also completely consistent with the recent analytic treatments of Brady and Hughes ||31[ 



and of Wiseman |32| , which show that there can be no decrease in the maximum equilibrium 
rest mass at first PN order, presuming tidal effects can be neglected. However, we cannot 
be sure that PN tidal effects would not change our results (see §|V1 C]). Nevertheless, our 



results are consistent with the relativistic integrations of Baumgarte et al. (1997), which 
automatically include all PN tidal terms. This agreement provides strong support of our 
analysis, at least in the case of corotating binaries. 
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C. Directions for Future Work 



Recently, Wilson, Mathews, and Marronetti [It] have presented numerical simulations 
which suggest that otherwise stable neutron stars may individually collapse to black holes 
when brought into a close binary orbit. Our results suggest that the maximum equilibrium 
mass increases as the orbit decays. In this section, however, we discuss two important 
caveats which leave open the possibility that the stars could collapse prior to the orbital 
plunge: (1) the PN tidal terms neglected in this paper may affect the maximum equilibrium 
mass, and (2) even if the maximum equilibrium mass does increase in binaries, some of the 
allowable equilibria may be unstable to collapse. 

The largest uncertainty in our results may come from neglected PN tidal terms. Chan- 
drasekhar and Chandrasekhar and Nutku |0| have derived the PN corrections to the 
energy and angular momentum for arbitrary fluid configurations, which could be used to 
determine the PN corrections for our binary star system. The tidal components of these 
corrections are absent from our present analysis. For example, consider the PN correction 
from Ref. @]'s equation (68) which is proportional to / d^x, where $ is the usual New- 
tonian gravitational potential. Inside a star, the potential $ can be decomposed into an 
internal contribution due to that star, and an external contribution ^^xt due to the 
orbital companion. The integral / p$|„j d^x then contributes to the self-gravitational PN 
energy Wpn, while the cross term / p^ext^intd^x gives a coupling as in Eqi- The integral 
/ p^l^f-d^x to lowest order helps give the point mass PN orbital energy like in Ref. @], and 
for finite-size stars this integral can be evaluated in a multipole expansion with the higher 
order terms accounting for PN tidal effects. A careful investigation of PN tidal effects seems 
worthwhile. 

We also note that we have not tested our equilibrium configurations for stability against 
collapse to black holes: even if the maximum equilibrium mass does increase as the orbit 
decays, this does not necessarily mean the maximum stable mass follows the same trend. The 
EFE energy variational approach provides a straightforward means of testing the stability 
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of equilibrium models, which could be exploited in future studies. In particular, the second 
order variation of the energy functional can be used to identify any instability to collapse. 
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APPENDIX A: POST-NEWTONIAN EXPANSION OF THE OV EQUATION 

In this appendix we expand the OV equation (see, eg. |^) to PN order and thereby derive 
the first PN corrections to the Newtonian Lane-Emden functions for an isolated spherical 
polytrope in equilibrium. Our results are useful for generating trial density functions in § |B 2| 
accurate to PN order. 

Let Poir) be the rest mass density of a spherically symmetric distribution of matter at 
the Schwarzschild radius r. Let m{r) be the enclosed total mass-energy, and define the 
quantity C, and the non-dimensional functions and z/(^) by 

Po = Poc^'^, r = a^, m = -ina^poci^^'^, (Al) 

where poc = Po(0) is the central rest mass density, so that -(9(0) = 1. In the Newtonian 
limit, the functions and z/(,^) reduce to the usual Lane-Emden functions 6 and 6'. Here 
a = R/^i is a scale factor, with R being the radius of the star in Schwarzschild coordinates 
and ^/ being defined by the condition = 0. For a polytropic equation of state 

p = Kpl^- (A2) 

the OV equations yield a hydrostatic equilibrium solution if (see, eg., Ref. fSSH ) 
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dd_ ^ t/(l + (n + l)g,^)(l-g,^"+^eH 

di (l + 2(n + l)gX) ' ^ ^ 

1 d{ve) 



-T{l + nqcd), (A4) 
,2 9c, (A5) 



oc 



C 

where we have defined 

Qc = 7 — ^YT- (A6) 
[n + l)c^ 

Evaluating the last equation of at the surface, we see that for an equilibrium configu- 
ration the scale factor a satisfies 

iira^Poc = -J^, (A7) 

where M is the total mass-energy and uj = In the Newtonian limit (c oo), q^. —>■ 0, 

d ^ 9 and u —>■ 6', where 6 is the usual Lane-Emden function. If we define the functions a 
and f3 as the PN components in the expansions 

^ = e + q,a + 0{q',), (A8) 
u = e' + q,P + 0{q''J, (A9) 

then equations (^) and (|A4|) become 



— = f3 + (n + 1)96' - r+^e - 2(n + (AlO) 



1 d{P^ 



2\ 



C2 AC =-nM (All) 

In deriving equations (|A10|) and (|A11|) we have used the Lane-Emden equation to cancel 



all terms of Newtonian order. We note that the functions a and (3 depend only on the 
polytropic index n and can be numerically integrated simultaneously with the Lane-Emden 
equation. The numerical integrations can be started slightly away from ^ = with the help 
of the approximations, for ^ << 1: 

a~-y(^ + 2), (A12) 
/? ^ (A13) 
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The surface occurs when = = 0(0) + Qc<^{^i) = (0 — + gc«(^i), to PN order, 
where we have Taylor expanded and used 9{^i) = to obtain the last equality. Therefore, 

e/ = ei-gc^ + 0(a, (A14) 

where ai = Similarly, Taylor expanding equation (|A9|) about the surface at ^ = 

gives 

1^1 = 0[ + + + 0{ql), (A15) 

where vj = and /?i = 

In the special case n = 1 it is well known that the solution to the Lane-Emden equation 
is ^ = sin(^)/^. Equations (|A10 ) and ( All ) can also be solved analytically when n = 1: 



a = 9'^ -9- 2^9' 9 - 3^^^^ - W ^9^ + 9 (^9' + ^"^9^ d^, (A16) 

P = 3^9' -9' + 299' - 3,^0^0' + 2^9'^ - W ^9' - + + ^2^3 .^17) 

Jo t, JO 

These equations provide a convenient check of our numerical integrations in Appendix B. 



APPENDIX B: POST-NEWTONIAN SELF-ENERGY TERMS 



We now derive the PN self energy terms of a spherically symmetric fluid which is instan- 



taneously at rest. Our approach generalizes the method of Ref. |^. Let us introduce the 
enclosed rest mass m' and the proper radial coordinate r', which are related to the enclosed 
total mass-energy m and the Schwarzschild radius r by 

dm' = poAnr'^ dr' = ( 1 - ) Po47rr^ dr = il - ( 1 + ^ j dm, (Bl) 

where u is the internal energy density. To PN order, the energy of such a system, excluding 



rest mass energy, is [see equations (6.9.3), (6.9.8), (6.9.10) and (6.9.18) of Ref. |^ 



pM j-Mo 

E = c^ dm-c^ dm' (B2) 
Jo Jo 
fMo fMo rn' 

= udm'~G —dm' + h + h + h + h + h, (B3) 
JO Jo r' 
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where Mo is the total rest mass, and 



Ii = — - / u—dm, (B4) 
c"' Jo r 

^2 M 



4=4r;^r„„„, (B6) 

JO r Jo 

h = ^ — -dm, B7 
Jo r Jo r 

h = 7 / — 3— / "^'^ u?^- (B8) 

C'^ Jo r** Jo 

Since we are working only to PN order, the quantities m and r can be replaced with m' and 
r' in Ji to Is, if desired. 

Note that equations ( [B2| ) and ( P3|) do not require that the fluid be in equilibrium: the 
density profile po and internal energy profile u need to be spherically symmetric, but are 
otherwise arbitrary. Let us now restrict our attention to fluids obeying a polytropic equation 
of state 

P = KpV^, u = nKpl''', (B9) 

where is a constant, but allow the density profile po to remain arbitrary for the moment. 
We wish to find the contribution to the energy beyond that given by the usual Newtonian 
expressions for the internal energy U and the self-gravitational energy W [see equations (|12D 
and dll)], ie. 

= E - (k.Kpli^Mo - k^Gpli^Ml'^) , (BIO) 

where it is understood that all terms of order P^N and higher are dropped. Inspection of 
equations (|B3|), ( P9| ) and ( |B10D reveals that /S.E is composed of two types of terms: those 



which depend on K and those which do not; we write 

AE = UpN + WpN, (Bll) 

where we have defined 
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u 



PN= u dm' - kiKpl'j'Mo + h + /3, (B12) 

JO 

WpN = ~G / —dm' + k2Gpli^M'J^ + l2 + h + h. (B13) 
JO r 

The quantity U + W + Upn + Wpn is the energy of a spherical symmetric configuration with 
an arbitrary density profile, to PN order. 

If we restrict our attention to a set of well-chosen density functions, then the integrals 
in equations ( [B12D and (|B13| ) can be evaluated. In particular, if we consider a family of 



density profiles of the form Po{r') oc f{r') [or equivalently Po{r) oc g{r)], where f{r') [or g{r)] 
is arbitrary, then 

UpN = -h%Kpt^'Mll\ (B14) 
Wp^ = -l2%p-i^Mll\ (B15) 

Here the dimensionless coefficients l\ and I2 are determined by evaluating the integrals in 
equations ([B12|) and ([B13D for the chosen family of density profiles. 



1. Lane-Emden Trial Density Functions 

As a concrete example, consider density profiles which are related to the non-dimensional 
Lane-Emden function Q by 

p'f{r')=p^r. r' = ^e, (B16) 

where i?' is the proper radius of the star and the central density poc does not necessarily 
have the value which would allow for equilibrium at Newtonian order. The superscript "LE" 
appearing in equation ( [B16| ) is a reminder that we are limiting ourselves to a single family 
of trial density profiles based upon Lane-Emden functions. Since density profiles of this 
form differ at least at PN order from the OV solution which minimizes the energy given by 



33 



equation (p^), there is no value of poc for which the fluid is in equilibrium to PN order. [| 
Using equation (|B16| ) in equations ( |B12| ) and ( |B13| ) we obtain simply 



= /^t'' 5/3 ' 



^2 - 2/3 7/3- 

~2h'oc J-yJ-o 



c 



From here it is straightforward to write and as integrals over Lane-Emden functions, 
and these coefficients therefore depend only on the polytropic index n. Equations ( pi7|) and 
( P18|) define the coefficients used to produce the short dashed curves in Figure |l|. 



In the following subsection we go through a considerable amount of work to derive /i 
and I2 for more accurate trial density functions. However, the coefficients given by equations 
( [B17|) and ( |B18| ) already represent a significant improvement over that used previously in the 



literature (compare the short-dashed and dot-dashed curves in Figure |I]). The reason for this 
improvement is that we have kept both of the PN terms [eqs. ( [B14[ ) and (|B15| )] which result 
naturally from the subtraction in equation ( [BlOl ). Previous analyses have used equilibrium 



relations to eliminate K in equation ( [B14D and forced UpN to scale like WpN, resulting in 



an energy functional which is valid only along equilibrium sequences and therefore should 
not be used in an energy variational treatment [see the discussion surrounding eq. ( P47|) 



below]. Nowhere in the derivation of equations ( P17|) and (P18|) have we assumed the star 
is in equilibrium. 

2. OV-Based Trial Density Functions 

The disadvantage of being restricted to Lane-Emden trial density functions is that the 
energy minimum of equation (|B3|) then differs at PN order from the actual energy minimum. 



We know that fluid conflgurations which obey the OV equation are precisely those necessary 



^The exception to this statement occurs for n = 0. In this case, both the purely Newtonian and 
the exact general relativistic density profile are of the form po =const. 
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to minimize the energy given by equation (B2). This suggests that we construct a family of 



trial density functions based upon the solution to the OV equation, namely 

Po(r) = Poc{e + ge«)", r = ^C, (B19) 

where R is the stellar radius in Schwarzschild coordinates, the central density poc is not 
necessarily the value required by the OV equation for equilibrium [ie. eqs. ( |A5D and (|A7|) 



need not be satisfied], and the function a is obtained by solving equations ( |A1CI|) and (|A11D . 
The expansion parameter qc is defined by equation (|A6|) with a = R/^i, ie. 

^|(n + l)c2 

An energy variational method based upon this family of trial density functions can identify 
the minimum of equation ( |B3D to PN order. There are an infinite number of other families 
which satisfy this criterion [for instance, we could generate one other such family by replacing 
Poc with its equilibrium value in eq. ( |B2CI| )]. However, the trial density functions defined by 
equations ( |B19D and ( |B2(J| ) are independent of the polytropic constant K and hence make 
it easy to track which terms do and do not depend on K. We note that the last equation of 
does not apply here, although by construction it is satisfied to PN order when poc has 
its equilibrium value. [To derive the enclosed total mass-energy even for non-equilibrium 
configurations, use m = J pAnr"^ dr with p = po{l + u/c^).] 

Our plan is now to evaluate li and I2 as defined by equations (|B14|) and ( |B15| ), using the 



family of trial density functions given in equation (pi9|) . We begin with the first integral 
appearing in equation ( |B12| ) and work to PN order: 

\(l^^.lL),,„ ,B22) 



\ rc^ (p- 

M / u 



u[l-—j dm- /i, (B23) 
where we have used equation (|Bl| ) and expanded. We now focus our attention on the integral 
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M / u\ rR 



V c^y "'0 



u[l- ^] dm = Ann J Pr"^ dr (B24) 
ni^pV"^^^ f\e + qcar'e dC (B25) 
nA>;/»l^ r (l + (n + 1)^1 e (B26) 



= nKpli-^^^ [^f^' e-^'e di+{n + f^' aO^e d^j ■ (B27) 

Note that since we are working only to PN order, it is justified to end all the integrals in 
this appendix at [where is defined by ^(^i) = 0]. For instance, since 6* ^ gc in the range 
6 < ^ < 6, and since - 6 ~ 

f^' r+^e' d^ = r+^e' d^ + = -^I'l^il + (^(g^')- (B28) 

JO Jo n 



Therefore equation ( |B27| ) becomes 



In order to write our final answer in the form of equation ( [Bll| ), we need to derive a 



relation between {AirR^poc/ ^f) and the total rest mass Mq. From our knowledge of Newtonian 
polytropes we expect {AtcR^Poc/Ci) — ^o/{^iWi\) + 0{qc), and we can obtain a relation 
accurate to PN order by evaluating Mq: 

rMo 

Mo = / dm' (B30) 
Jo 



, ,„+g,ar + (B32) 



( lo' + (nar-^e' - (n + 1)^"^'^') d^^ . (B34) 



Making use of 



/ ' di = ei\e[\ + o(gr^), / ' e-e'e d^ = - — e^w + o(gr^), (b35) 

JO JO b — n 
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in equation (P34|) gives 



Mo 
Mo 



1 



nqc 

■:'2\ai 



ei\e'i\ Jo 



5 — n 



^ + ^^Ton + E Z 



5 — n 



(B36) 
(B37) 



where /3i = /3(G)- We evaluated the integral in equation (p336|) with the help of equations 
( IATiI) and ( [B28|) . Equation ( p37|) is our desired relation between {AirR^poc/^i) and Mo. 
Using equations (|B37| ) and (|B29|) in equation ( P23|) gives 

(n-3)(n + 1) 



5 — n 



(B38) 



We note that the lowest order term is the usual internal energy contribution for a pure 
polytrope, kiK pUj" Mo, as expected; the first order corrections will help give us our expression 
for the coefficient h. Using equation ( |B38| ) and the definition of ki = n{n + l)G|^^'J/(5 — n) 
in equation (P12|) gives 



Kpli^Monin + / /^iCi^l^^ 



(5-n)2 / 

or, after using equations ( [B20| ) and ( |B37| ) to lowest order and simplifying, 

(47r)V3 



(B39) 



(n -3)(n + 1) 



(5 — n) 



(5 - ny 



(B40) 



The relation^ 



^To prove this relation, begin by non-dimensionalizing the double integral I3 and writing it in 
terms of the single integrals Ii and I2 in non-dimensional form, following the strategy of Exercise 
6.19 in Ref. |^0|. Then with repeated integration by parts and use of the Lane-Emden equation, 
the two remaining integrals can be written in terms of each other, and the identity is obtained. 
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d^e-U'0'''-'ed^ = ^^M'-^ Tee'e^'^'d^ (B4i) 

JO 5 — n 6 Jo 

can be used to reduce the double integral in equation (P4CI|) to a single integral: 

(4^)V3 J „-l /•«.,3^,_, 



If desired, the identity [which can be proved from equations ( |A1CI| ) and ( |A11| ) and appropriate 
integrations by parts] 



Jo 1 — n 6 ^ ^ Jo 



(B43) 



could be used to further manipulate equation (P42|) when n ^ 1, while for n = 1 this integral 
can be evaluated using equation ([A16|) . However, we find it more convenient to work with 
equation (|B42| ) directly. 

Equation ( P42| ) is our final expression for the coefficient li, and was obtained by evalu- 
ating equation (pi2|) with the OV-based trial density functions defined in equation ( |B19| ). 
We could proceed to calculate I2 from equation ( |B13| ) in a similar manner. However, it is 
more straightforward to take the following approach. 

Rearranging equations ( piO| ) and ( pil| ) yields the PN accurate energy functional for an 
isolated spherical star, not necessarily in equilibrium, obeying a polytropic equation of state: 

E = hKpli^M, - k^Gpl'^'Ml'^ - h%Kpt^'Ml'^ - l^^p^l^Ml'^ (B44) 

c c 

Along an equilibrium sequence we have = dE/dpoc, which gives, to lowest order, 

Kpl!:^ = "^GMTpI/I (B45) 

Indeed, equation (|B45|) is the same relation obtained by writing the Newtonian equation of 
hydrostatic equilibrium in terms of Lane-Emden functions and simplifying. To PN order, 
the energy along an equilibrium sequence therefore satisfies 



2 



E^^ = hKpl'-M^ - hGpli'Ml'-^ - (/i^ + k)^pl/'Ml'\ (B46) 
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Since here we are focusing on the equihbrium sequence, the first variational derivative of 
E with respect to the density functional must vanish; therefore, any variation in the density 
profile of PN order causes a variation in the energy of order P^N, so that it is sufficient to use 
Lane-Emden profiles when evaluating the Newtonian order integrals in AE. We therefore 
find that along an equilibrium sequence AE = /1 + /2 + /3 + /4 + /5. The Newtonian equation 
of hydrostatic equilibrium, equation (|B45|) , can then be used to evaluate Ji through J5 in 



this case. These integrals have been carried out by Ref. (see their §6.9), so that 

E^'^ = hKpli-M, - hGpli'M'J' - k^pli^Ml/^ (B47) 



c2 



where (see eq. [6.9.31] of 







^ ^ - ^ - ^^^±^1^ eO'r-^ d^] . (B49) 



M\6[\Y/^ \2{5-ny^' ^' 6(n+ 1) Jo 

By integrating the OV equation exactly and extracting coefficients from the resulting equi- 
librium sequence energy E^'^{poc), we have numerically verified that equation ( |B47| ) holds 
to PN order. We emphasize that equation ( |B46| ) or ( |B47| ), in contrast to equation ( |B44D , 
should not be used as energy functionals in an energy variational treatment, since they are 
valid only for equilibrium configurations^; this has not been noticed in previous analyses. 
Comparing equations (|B46D and (|B47| ) we see immediately 



l2 = k-h^ (B50) 



3 



In the special case n = 0, the self energy correction Upn = 0, and the use of equation (B47) as 



an energy functional does yield correct answers to first PN order. Equation ( |B47] ) also yields PN 

2 /3 

accurate answers when n = 3; this is because Upn and Wpn both then scale as poc and because 
the equilibrium relation between K and Mq is independent of poc at lowest order [see eq. (B45)]. 



However, for n ^ and n 7^ 3, use of equation ( p47 ) as an energy functional gives errors at first 
PN order. 
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+ 



- lOn - 5 



2(5 -n)2 




(B51) 



(B52) 



Equations (P42| ) and ( P52| ) are our final expressions for tlie coefficients li and I2, and were 



used to generate the long dashed curves in Figure |I|. The Newtonian curves in this figure, as 



M which is not accompanied by a maximum in the rest mass Mo] this is because along these 
sequences dM/dpoc = {dM/ dMo){dMo/dpoc) vanishes due to dM/dMo = 0. The remaining 
sequences have maxima in M and Mo at the same poc, where = dM/dpoc = dMo/dpoc, and 
these turning points mark the onset of radial instability. 

The derivation of the coefficients li and I2 is formally invalid if n = 0, since then the 
first relation in equation (P35| ) is accurate only to Newtonian order. However, for n = one 
can analytically derive the coefficients exactly to be Zi = and I2 = 3(47r/3)^/^/70, and the 
numerically obtained values given by equations (P42|) and ( |B52|) do indeed approach these 
limiting values as n approaches zero. See Table |I| for numerical values of these coefficients 
for various n. We note that the coefficient li is fairly well approximated simply by li ~ n. 
Also note that for n > 0.5, li is roughly 10 or more times larger in magnitude than I2, so 
that we expect Upn to be significantly larger than Wpn in magnitude for all but the stiffest 
equations of state. 



well as the n = 1 curves which use equation ( |B47| ), have a maximum in the gravitational mass 
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TABLES 

TABLE L The polytropic coefficients to tlie PN energy corrections of an isolated spherical star, 
as determined from equations ( B42 ) and ( B52 ). 



n h h 

0.0 0.000000 0.111365 

0.1 0.094949 0.097502 

0.2 0.191634 0.084840 

0.3 0.289737 0.073070 

0.4 0.388996 0.061974 

0.5 0.489198 0.051391 

0.6 0.590166 0.041201 

0.7 0.691752 0.031312 

0.8 0.793835 0.021649 

0.9 0.896311 0.012154 

1.0 0.999096 0.002779 

1.1 1.102118 -0.006516 

1.2 1.205321 -0.015767 

1.3 1.308657 -0.025003 

1.4 1.412088 -0.034252 

1.5 1.515586 -0.043538 
2.0 2.033637 -0.091306 
2.5 2.553451 -0.143291 
3.0 3.080363 -0.202770 
3.5 3.628218 -0.275041 
4.0 4.231742 -0.371100 
4.5 5.003087 -0.523069 
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TABLE II. Isolated spherical polytropic models, as calculated from the PN-accurate energy 
functional ( B44| ). Here n is the polytropic index; KpJ^ jc? is the central density parameter; M is 



the gravitational mass; K is the stellar radius in Schwarzschild coordinates as calculated from eq. 
dll); Mq is the rest mass; and Kq is the Newtonian radius defined by eq. (|4^). The value of Kpoc/f? 
marked by an asterisk (*) identifies the star with maximum rest mass Mq and gravitational mass 



M for that value of n. 
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TABLE III. Irrotational Equilibrium Sequences. The dimensionless quantities fi, J, and E are 



defined in equation (43). 
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8814 


0.373(-3) 


0.3084 


1.2603 


-1.2698 


1.0014 


2.6 


2.920 


oo 


0.8322 





8537 


0.724{-3) 


0.3340 


1.2401 


-1.2754 


1.0023 


2.344* 


2.759 


oo 


0.7716 





8049 


0.178(-2) 


0.3681 


1.2283 


-1.2790 


1.0042 


2.0 


2.623 


oo 


0.6498 





7045 


0.635(-2) 


0.4090 


1.2651 


-1.2663 


1.0107 










GM/Rc^ 


= 0.2, n = 0.5 










6.0 


5.173 


13.36 


0.9845 





9848 


0.424(-6) 


0.1146 


1.5553 


-1.2860 


0.8533 


5.0 


4.344 


11.22 


0.9732 





9737 


0.210(-5) 


0.1458 


1.4727 


-1.2953 


0.8534 


4.0 


3.537 


9.14 


0.9475 





9495 


0.146(-4) 


0.1930 


1.3956 


-1.3064 


0.8536 


3.5 


3.150 


8.14 


0.9221 





9259 


0.457(-4) 


0.2260 


1.3647 


-1.3120 


0.8539 


3.2 


2.928 


7.56 


0.8989 





9049 


0.972{-4) 


0.2499 


1.3520 


-1.3146 


0.8542 


3.0 


2.787 


7.20 


0.8784 





8866 


0.166(-3) 


0.2679 


1.3477 


-1.3155 


0.8547 


2.922* 


2.734 


7.06 


0.8690 





8782 


0.207(-3) 


0.2753 


1.3473 


-1.3156 


0.8549 


2.8 


2.652 


6.85 


0.8524 





8636 


0.293{-3) 


0.2875 


1.3485 


-1.3153 


0.8553 


2.6 


2.527 


6.53 


0.8192 





8345 


0.534(-3) 


0.3090 


1.3566 


-1.3132 


0.8563 


2.0 


2.258 


5.83 


0.6525 





6899 


0.396(-2) 


0.3794 


1.4748 


-1.2770 


0.8653 


1.736 


2.245 


5.80 


0.5351 





5835 


0.103(-1) 


0.4036 


1.6258 


-1.2255 


0.8772 










GM/Rc^ 


= 0.25, n = 0.5 










6.0 


4.878 


10.32 


0.9836 





9838 


0.437(-6) 


0.1171 


1.4952 


-1.3223 


0.8042 


5.0 


4.099 


8.67 


0.9716 





9721 


0.213(-5) 


0.1479 


1.4296 


-1.3297 


0.8043 


4.0 


3.340 


7.07 


0.9446 





9465 


0.144{-4) 


0.1938 


1.3749 


-1.3377 


0.8045 


3.5 


2.978 


6.30 


0.9182 





9218 


0.444(-4) 


0.2256 


1.3588 


-1.3406 


0.8048 


3.280* 


2.825 


5.98 


0.9014 





9064 


0.760{-4) 


0.2421 


1.3566 


-1.3410 


0.8051 


3.2 


2.770 


5.86 


0.8943 





8998 


0.931(-4) 


0.2486 


1.3570 


-1.3409 


0.8053 


3.0 


2.638 


5.58 


0.8733 





8808 


0.157(-3) 


0.2658 


1.3612 


-1.3400 


0.8057 


2.8 


2.512 


5.31 


0.8470 





8571 


0.274(-3) 


0.2847 


1.3717 


-1.3375 


0.8065 


2.6 


2.396 


5.07 


0.8136 





8274 


0.493{-3) 


0.3053 


1.3913 


-1.3324 


0.8076 


2.0 


2.145 


4.54 


0.6490 





6822 


0.352(-2) 


0.3748 


1.5575 


-1.2825 


0.8174 


1.682 


2.146 


4.54 


0.5069 





5517 


0.112{-1) 


0.4034 


1.7850 


-1.2052 


0.8341 
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GMjRi? = 0.2, n = 1.0 



6.0 


4.026 


13.72 


0.9894 


0.9895 


0.144 


-6) 


0.1713 


1.4318 


-1.3369 


0.6663 


5.0 


3.373 


11.50 


0.9816 


0.9819 


0.719 


-6) 


0.2188 


1.3528 


-1.3502 


0.6663 


4.0 


2.732 


9.31 


0.9638 


0.9648 


0.506 


-5) 


0.2917 


1.2766 


-1.3668 


0.6667 


3.5 


2.422 


8.25 


0.9459 


0.9480 


0.161 


-4) 


0.3434 


1.2431 


-1.3759 


0.6672 


3.2 


2.241 


7.64 


0.9294 


0.9327 


0.346 


-4) 


0.3813 


1.2266 


-1.3811 


0.6679 


3.0 


2.125 


7.24 


0.9147 


0.9193 


0.598 


-4) 


0.4098 


1.2182 


-1.3840 


0.6686 


2.8 


2.013 


6.86 


0.8957 


0.9022 


0.107 


-3) 


0.4413 


1.2127 


-1.3860 


0.6698 


2.660* 


1.939 


6.61 


0.8792 


0.8874 


0.163 


-3) 


0.4651 


1.2115 


-1.3865 


0.6710 


2.6 


1.908 


6.50 


0.8712 


0.8803 


0.197 


-3) 


0.4757 


1.2117 


-1.3864 


0.6717 


2.0 


1.664 


5.67 


0.7405 


0.7667 


0.159 


-2) 


0.5884 


1.2636 


-1.3615 


0.6889 


1.721 


1.636 


5.58 


0.6317 


0.6720 


0.481 


-2) 


0.6266 


1.3521 


-1.3146 


0.7143 










GM/Rc'^ 


= 0.25, n = 


= 1.0 










6.0 


3.398 


10.87 


0.9890 


0.9891 


0.149 


-6) 


0.2128 


1.3490 


-1.3952 


0.5621 


5.0 


2.850 


9.11 


0.9808 


0.9811 


0.739 


-6) 


0.2700 


1.2860 


-1.4083 


0.5626 


4.0 


2.316 


7.41 


0.9623 


0.9633 


0.512 


-5) 


0.3557 


1.2303 


-1.4232 


0.5647 


3.5 


2.063 


6.60 


0.9438 


0.9458 


0.160 


-4) 


0.4144 


1.2099 


-1.4299 


0.5675 


3.2 


1.919 


6.14 


0.9269 


0.9301 


0.342 


-4) 


0.4559 


1.2024 


-1.4327 


0.5706 


3.0 


1.827 


5.84 


0.9118 


0.9162 


0.587 


-4) 


0.4863 


1.2006 


-1.4334 


0.5737 


2.987* 


1.822 


5.83 


0.9107 


0.9152 


0.608 


-4) 


0.4884 


1.2006 


-1.4335 


0.5739 


2.8 


1.741 


5.57 


0.8926 


0.8987 


0.104 


-3) 


0.5189 


1.2023 


-1.4327 


0.5779 


2.6 


1.662 


5.32 


0.8678 


0.8765 


0.190 


-3) 


0.5531 


1.2089 


-1.4296 


0.5836 


2.0 


1.500 


4.80 


0.7382 


0.7628 


0.149 


-2) 


0.6538 


1.2835 


-1.3897 


0.6194 


1.688 


1.514 


4.84 


0.6158 


0.6556 


0.515 


-2) 


0.6787 


1.3948 


-1.3250 


0.6631 
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TABLE IV. Irrotational Equilibrium Sequences: Energy Terms 





11 4- IT' 


W + W 




vvt 




_ 








E 














Kjr IVl 1 JLt- 


— 0.2, n = 


1.0 
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n nni 

-U.UUIO 


n AQV/I 
U.Uo 1 4 


-i.ooDy 


o.u 


1 .014Z 


-2, 


1701 


U.UUUU 


fi nnni 

-U.UUUl 


n 1 zi An 

-U. I'-l^iU 


n nni 9 


n '^Q9A 


n nni 

-U.UUIO 


n nzif^n 

U.U40U 


- 1 . uz 




1 5121 


-2, 


1689 


n nnnn 

U.UUUU 


n nnn"^ 


n 1 71 4 


n nni 


n i^Qi 
-vj.oy lo 


n nni 

-U.UU-LO 


U.UOUZ 


- -L.OUUo 


O.O 


1 .OUoD 


-2, 


,1669 


A nnnn 

U.UUUU 


A nni n 
-u.uuiu 


n 1 R7Q 

-u. lo 1 y 


n nm fi 

-U.UUIO 


-u.ooyo 


n nni a 

-U.UUIO 


n nfi'^s 

U.UOoo 


1 '37'^Q 

-i.o / oy 




1.0U4Z 


-2. 


,1644 


n nnni 

U.UUUl 


n nni fi 
-u.uuio 


n 1 QS1 
-u. lycsi 


n nni 7 


n "=1879 

-u.oo / z 


n nm 4 

-U.UUl'i 


n nfiQ9 


1 '^81 1 


3.0 


1.4993 


-2. 


,1616 


0.0001 


-0.0023 


-0.2047 


-0.0018 


-0.5847 


-0.0014 


0.0731 


-1.3840 


2.8 


1.4917 


-2. 


,1572 


0.0002 


-0.0034 


-0.2106 


-0.0018 


-0.5807 


-0.0014 


0.0774 


-1.3860 


2.660* 


1.4839 


-2. 


,1526 


0.0004 


-0.0045 


-0.2140 


-0.0019 


-0.5767 


-0.0014 


0.0804 


-1.3865 


2.6 


1.4797 


-2, 


1502 


0.0004 


-0.0051 


-0.2151 


-0.0019 


-0.5745 


-0.0014 


0.0817 


-1.3864 


2.0 


1.3752 


-2, 


,0864 


0.0033 


-0.0182 


-0.2027 


-0.0027 


-0.5211 


-0.0014 


0.0925 


-1.3615 


1.696 


1.2219 


-1, 


,9844 


0.0106 


-0.0349 


-0.1615 


-0.0041 


-0.4451 


-0.0013 


0.0906 


-1.3082 





TABLE V. ImAi 


iliovui] E('|uilil)rium S 


H'f|U(nic('s: Ang 


:ular MouKuilum Terms 




r/a\ 


Jh 


Js + J's 


JpN 


Joi + J'oj 


J 






GMIBx? = 


0.2, n = 1.0 






6.0 


1.5437 


0.0000 


-0.0017 


-0.1102 


1.4318 


5.0 


1.4524 


0.0000 


-0.0008 


-0.0988 


1.3528 


4.0 


1.3623 


0.0001 


0.0006 


-0.0864 


1.2766 


3.5 


1.3213 


0.0002 


0.0015 


-0.0799 


1.2431 


3.2 


1.3000 


0.0005 


0.0022 


-0.0760 


1.2266 


3.0 


1.2882 


0.0007 


0.0026 


-0.0734 


1.2182 


2.8 


1.2794 


0.0012 


0.0030 


-0.0709 


1.2127 


2.660* 


1.2757 


0.0017 


0.0032 


-0.0692 


1.2115 


2.6 


1.2749 


0.0021 


0.0033 


-0.0685 


1.2117 


2.0 


1.3112 


0.0130 


0.0032 


-0.0638 


1.2636 


1.696 


1.3890 


0.0388 


0.0008 


-0.0647 


1.3638 
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TABLE VI. Equilibrium Corotational Sequences. The dimensionless quantities i7, J, E, and C 



are defined in equation ([43|). 



r/ai 


r/Ro 


rc^/Gm 


0.2/0.1 


03/01 


Ts/\W\ 


n 




J 




E 


R/ Ro 


C 










GM/R 


= (Newtonian treatment), n 


= 0.5 










6.0 


6.073 


oo 





9874 


0.9794 


0.221(-2) 


0.1092 


1 


8057 


-1 


1904 


1.0009 


-0.0623 


5.0 


5.105 


oo 





9784 


0.9651 


0.376(-2) 


0.1418 


1 


6809 


-1 


2039 


1.0016 


-0.0816 


4.0 


4.162 


oo 





9585 


0.9348 


0.711(-2) 


0.1929 


1 


5603 


-1 


2211 


1.0031 


-0.1130 


3.5 


3.710 


oo 





9390 


0.9065 


0.103(-1) 


0.2299 


1 


5072 


-1 


2307 


1.0046 


-0.1369 


3.2 


3.450 


oo 





9214 


0.8819 


0.131(-1) 


0.2571 


1 


4810 


-1 


2363 


1.0061 


-0.1554 


3.0 


3.283 


oo 





9058 


0.8609 


0.156(-1) 


0.2777 


1 


4674 


-1 


2394 


1.0074 


-0.1701 


2.8 


3.124 


oo 





8860 


0.8352 


0.186(-1) 


0.3005 


1 


4583 


-1 


2417 


1.0092 


-0.1871 


2.625* 


2.992 


oo 





8641 


0.8079 


0.219(-1) 


0.3221 


1 


4553 


-1 


2425 


1.0112 


-0.2042 


2.6 


2.974 


oo 





8605 


0.8035 


0.224(-l) 


0.3254 


1 


4554 


-1 


2424 


1.0115 


-0.2068 


2.115** 


2.677 


oo 





7609 


0.6909 


0.367(-l) 


0.3915 


1 


4949 


-1 


2298 


1.0218 


-0.2696 


2.0 


2.630 


oo 





7256 


0.6543 


0.416(-1) 


0.4066 


1 


5208 


-1 


2208 


1.0259 


-0.2880 












GM/R 


= 0.2, n 


= 0.5 














6.0 


5.189 


13.40 





9848 


0.9782 


0.175(-2) 


0.1141 


1 


6096 


-1 


2832 


0.8541 


-0.0524 


5.0 


4.365 


11.27 





9739 


0.9634 


0.284(-2) 


0.1449 


1 


5423 


-1 


2907 


0.8546 


-0.0671 


4.0 


3.565 


9.21 





9499 


0.9324 


0.503(-2) 


0.1911 


1 


4893 


-1 


2983 


0.8558 


-0.0899 


3.5 


3.182 


8.22 





9269 


0.9039 


0.699(-2) 


0.2232 


1 


4756 


-1 


3007 


0.8571 


-0.1067 


3.345* 


3.067 


7.92 





9172 


0.8920 


0.780(-2) 


0.2348 


1 


4746 


-1 


3009 


0.8576 


-0.1129 


3.2 


2.962 


7.65 





9065 


0.8792 


0.867(-2) 


0.2465 


1 


4756 


-1 


3007 


0.8582 


-0.1194 


3.0 


2.821 


7.29 





8887 


0.8584 


O.lOl(-l) 


0.2640 


1 


4809 


-1 


2995 


0.8593 


-0.1294 


2.842** 


2.714 


7.01 





8716 


0.8387 


0.115(-1) 


0.2790 


1 


4893 


-1 


2975 


0.8604 


-0.1383 


2.8 


2.687 


6.94 





8666 


0.8329 


0.119(-1) 


0.2832 


1 


4922 


-1 


2968 


0.8608 


-0.1408 


2.6 


2.560 


6.61 





8388 


0.8017 


0.141(-1) 


0.3042 


1 


5117 


-1 


2918 


0.8627 


-0.1542 


1.734 


2.231 


5.76 





5976 


0.5539 


0.357(-l) 


0.4054 


1 


8183 


-1 


1935 


0.8899 


-0.2491 












GM/Rc^ = 0.25, n 


= 0.5 














6.0 


4.892 


10.35 





9839 


0.9778 


0.160(-2) 


0.1167 


1 


5459 


-1 


3196 


0.8050 


-0.0492 


5.0 


4.117 


8.71 





9723 


0.9627 


0.257(-2) 


0.1470 


1 


4941 


-1 


3255 


0.8055 


-0.0625 


4.0 


3.365 


7.12 





9470 


0.9312 


0.447(-2) 


0.1921 


1 


4606 


-1 


3303 


0.8066 


-0.0830 


3.706* 


3.151 


6.67 





9341 


0.9157 


0.537(-2) 


0.2096 


1 


4580 


-1 


3307 


0.8073 


-0.0912 


3.5 


3.006 


6.36 





9227 


0.9022 


0.615(-2) 


0.2232 


1 


4596 


-1 


3304 


0.8078 


-0.0979 


3.230** 


2.820 


5.96 





9038 


0.8801 


0.741(-2) 


0.2432 


1 


4674 


-1 


3288 


0.8089 


-0.1079 


3.2 


2.800 


5.92 





9013 


0.8773 


0.757(-2) 


0.2456 


1 


4688 


-1 


3285 


0.8090 


-0.1091 


3.0 


2.668 


5.64 





8828 


0.8562 


0.878(-2) 


0.2624 


1 


4813 


-1 


3258 


0.8101 


-0.1180 


2.8 


2.542 


5.38 





8598 


0.8305 


0.103(-1) 


0.2809 


1 


5010 


-1 


3211 


0.8116 


-0.1282 


2.6 


2.425 


5.13 





8312 


0.7990 


0.122(-1) 


0.3012 


1 


5304 


-1 


3137 


0.8136 


-0.1401 


1.679 


2.131 


4.51 





5634 


0.5264 


0.338(-l) 


0.4053 


1 


9611 


-1 


1760 


0.8463 


-0.2354 
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GM/Rx? = 0.2, n = 1.0 



6.0 


4.045 


13.79 


0.9895 


0.9847 


0.127 


-2) 


0.1702 


1 


4714 


-1 


3338 


0.6683 


-0.0373 


5.0 


3.398 


11.58 


0.9820 


0.9741 


0.208 


-2) 


0.2166 


1 


4040 


-1 


3450 


0.6697 


-0.0479 


4.0 


2.769 


9.44 


0.9651 


0.9516 


0.373 


-2) 


0.2867 


1 


3461 


-1 


3574 


0.6728 


-0.0645 


3.5 


2.466 


8.41 


0.9486 


0.9303 


0.523 


-2) 


0.3354 


1 


3260 


-1 


3628 


0.6759 


-0.0769 


3.2 


2.292 


7.81 


0.9337 


0.9117 


0.652 


-2) 


0.3706 


1 


3195 


-1 


3647 


0.6787 


-0.0863 


3.056* 


2.211 


7.53 


0.9246 


0.9005 


0.729 


-2) 


0.3893 


1 


3186 


-1 


3650 


0.6805 


-0.0916 


3.0 


2.180 


7.43 


0.9206 


0.8957 


0.762 


-2) 


0.3969 


1 


3188 


-1 


3650 


0.6813 


-0.0938 


2.8 


2.073 


7.06 


0.9041 


0.8759 


0.899 


-2) 


0.4255 


1 


3220 


-1 


3638 


0.6848 


-0.1025 


2.6 


1.972 


6.72 


0.8831 


0.8512 


0.107 


-1) 


0.4565 


1 


3307 


-1 


3605 


0.6895 


-0.1126 


2.570** 


1.957 


6.67 


0.8795 


0.8470 


0.110 


-1) 


0.4613 


1 


3326 


-1 


3597 


0.6903 


-0.1143 


1.730 


1.703 


5.81 


0.6886 


0.6420 


0.276 


-1) 


0.5933 


1 


5182 


-1 


2710 


0.7500 


-0.1913 










GM/Bx 


.2=0. 


25, n 


= 1.0 














6.0 


3.460 


11.06 


0.9891 


0.9845 


0.122 


-2) 


0.2077 


1 


3903 


-1 


3908 


0.5715 


-0.0339 


5.0 


2.926 


9.36 


0.9813 


0.9738 


0.197 


-2) 


0.2607 


1 


3385 


-1 


4012 


0.5763 


-0.0435 


4.0 


2.409 


7.70 


0.9638 


0.9510 


0.350 


-2) 


0.3381 


1 


2990 


-1 


4113 


0.5849 


-0.0585 


3.5 


2.163 


6.92 


0.9468 


0.9296 


0.487 


-2) 


0.3901 


1 


2895 


-1 


4142 


0.5921 


-0.0698 


3.371* 


2.101 


6.72 


0.9408 


0.9222 


0.534 


-2) 


0.4054 


1 


2890 


-1 


4144 


0.5945 


-0.0733 


3.2 


2.022 


6.47 


0.9315 


0.9108 


0.606 


-2) 


0.4267 


1 


2900 


-1 


4140 


0.5982 


-0.0784 


3.0 


1.933 


6.18 


0.9181 


0.8947 


0.707 


-2) 


0.4534 


1 


2942 


-1 


4124 


0.6033 


-0.0853 


2.839** 


1.864 


5.96 


0.9049 


0.8791 


0.805 


-2) 


0.4763 


1 


3006 


-1 


4098 


0.6083 


-0.0917 


2.8 


1.848 


5.91 


0.9013 


0.8749 


0.832 


-2) 


0.4820 


1 


3026 


-1 


4090 


0.6097 


-0.0934 


2.6 


1.769 


5.66 


0.8799 


0.8502 


0.990 


-2) 


0.5123 


1 


3167 


-1 


4029 


0.6176 


-0.1029 


1.700 


1.596 


5.10 


0.6739 


0.6301 


0.269 


-1) 


0.6335 


1 


5440 


-1 


2852 


0.7053 


-0.1832 
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TABLE VII. Physical quantities at the ISCO for selected irrotational sequences. Here 
Mi,4 = M/{1AMq) and Mf|^ = M'"^^/(1.6Mq), where M'"^^ is the maximum value of the 
isolated gravitational mass M. 



Ref. [n|'s Results 



MiAfisco M^l^'fisco 



Ml. 4f ISCO 



n 


GM/Rc^ 


GM/Roc^ 


JcjGMl 


Mil Ml 




aa/ai 


[Hz] 


[Hz] 


02/01 




[Hz] 


0.5 


0.125 


0.1145 


3.7178 


3.387 


0.821 


0.839 


258 


596 








0.5 


0.1380 


0.125 


3.5656 


3.326 


0.828 


0.845 


289 


597 


0.830 


0.850 


279 


0.5 


0.1941 


0.1667 


3.0998 


3.072 


0.865 


0.875 


404 


576 


0.857 


0.871 


399 


0.5 


0.2 


0.1707 


3.0622 


3.046 


0.869 


0.878 


413 


572 








0.5 


0.2481 


0.2 


2.7997 


2.841 


0.900 


0.905 


471 


531 


0.880 


0.891 


488 


0.5 


0.25 


0.2010 


2.7904 


2.833 


0.901 


0.906 


473 


529 








1.0 


0.125 


0.1008 


3.6758 


3.463 


0.829 


0.846 


297 


362 








1.0 


0.2 


0.1332 


3.1631 


3.228 


0.878 


0.887 


474 


438 








1.0 


0.25 


0.1405 


3.0364 


3.154 


0.911 


0.915 


542 


475 









TABLE VIII. Orbital frequency at the ISCO for selected corotating sequences. 
GMjRc^ M^.^ffico [Hz] Mi.4/,t^o [Hz] MJ-IV/Tco \M Mf^fV/lco [Hz] 



0.5 0.125 220 267 508 617 

0.5 0.2 353 419 488 580 

0.5 0.25 409 475 458 531 

1.0 0.125 248 300 302 366 

1.0 0.2 397 471 367 434 

1.0 0.25 450 529 394 463 
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FIGURES 

FIG. 1. We plot (a) the rest mass Mq and (b) the gravitational mass M = Mq + E/c? for 
isolated spherical stars as a function of the central density parameter Kp^J^ /(?. Here mass is 
in units c3-n(^-3/2^n/2 rpj^g solid curve was obtained by integrating the OV equation without 
approximation. The remaining curves were found by solving the equilibrium equation = dE/dpoc 
for various energy functionals E. The dotted-and-dashed curve results from the use of equation 
( p47D as the energy functional. The remaining three curves use different values of the coefficients 



li and I2, via equation ( B44| ): The dotted curve is the purely Newtonian result {li = ^2 = 0); the 



short-dashed curve uses the coefficients from equations ( B17 ) and ( B18 ); the long-dashed curve 



uses our standard coefficients from equations ( B42 ) and ( |B52 ) 



FIG. 2. The total energy E^ relative to its value E^ at infinite separation, as a function 
of the orbital frequency / = J7/27r for selected n = 0.5 irrotational sequences. The thick solid 
curves represent our FN sequences with various GM/Rc^, where M and R are the isolated neutron 
star gravitational mass and radius in Schwarzschild coordinates, respectively. The dashed curves 
represent purely Newtonian results, with the thin curve corresponding to the point mass sequence 
and the thick curve corresponding to the EFE treatment. The minima of these curves mark the 
ISCO, inside of which the orbit is dynamically unstable. The sequences terminate when the stars 
touch. 



FIG. 3. Same as Fig. ^, but for an n = 1 irrotational sequence. 

FIG. 4. The dependence of GM/Rc^ (as calculated in isolation) on the orbital frequency fisco 
at the ISCO for n = 0.5 and n = 1 irrotational sequences. For each n, the curves are terminated 
at the value GM/Rc^ corresponding to the maximum mass. 

FIG. 5. Same as Fig. |2|, but for an n = 0.5 corotating sequence. The minima of these curves 
mark the onset of the secular instability in the orbit. 
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FIG. 6. Same as Fig. ||, but for an n = 1 corotating sequence. 

FIG. 7. The rest mass Mq, normalized to the maximum rest mass of a star at infinite separation 
j^oo,max^ as a function of the central density parameter Qc = Kpic^ /c^ for an n = 1 corotating 
binary system and at three different values of the orbital angular frequency /: / = (solid curve), 
Mmaxj ^ (long dashed curve), and Mfg"""/ = 463 Hz (short dashed curve). The dotted 

continuations of the dashed curves represent those binaries inside the innermost dynamically stable 
orbit. 

FIG. 8. The maximum equilibrium rest mass M^^^ as a function of the orbital angular 
frequency / for corotating (dashed curves) and irrotational (solid curves) sequences. The curves 
are terminated at the values of / for which stars of rest mass Mq = M^'™^^ acquire a dynamically 
unstable orbit. 

FIG. 9. As in Fig. ^, we plot the binding energy versus orbital frequency / for various 
n = 1 corotating sequences: logq^c = —1-15 (solid curve), —0.87 (dashed curve), and —0.59 (dotted 
curve). The data points represented by crosses are from the computations of Baumgarte et al., for 
the same values of Qc- 
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